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A SEMI-DIRECT SOLVER FOR COMPRESSIBLE THREE-DIMENSIONAL ROTATIONAL FLOW 


Sin-Chung Chang and John J. Adamczyk 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 


SUMMARY 

An iterative procedure is presented for solving steady inviscid 3-D sub- 
sonic rotational flow problems. The procedure combines concepts from classical 
secondary flow theory with an extension to 3-D of a novel semi-direct Cauchy- 
Riemann solver. It is developed for generalized coordinates and can be exer- 
cised using standard finite difference procedures. The stability criterion of 
the iterative procedure is discussed along with its ability to capture the evo- 
lution of inviscid secondary flow in a turning channel. 


INTRODUCTION 

Quite often in fluid flows, one observes a complex cyclic flow pattern on 
planes normal to the primary flow direction. This flow structure is often 
referred to as "secondary flow." It is attributed to the streamwise component 
of vorticity that is generated by the deflection of the primary shear flow due 
to a boundary surface. Thus the primary driving force behind this phenomenon 
is the centrifugal pressure gradient in the main flow. The role of viscosity 
is of secondary importance. Practical examples of hardware in which one f.inds ? 
such flows are turbomachinery blade rows, transition diffuser ducting and fan- 
jet exhaust mixers. It is well known that the prediction and control of sec- 
ondary flow within these devices can lead to a substantial increase in their 
aerodynamic performance. 

One of the first theoretical efforts at describing the generation of sec- 
ondary flow was by Squire and Winter (ref. 1). This paper described the gen- 
eration of secondary flows in a turning channel, with the fluid assumed to be 
inviscid. Later Hawthorne (refs. 2 and 3) extended the analysis of Squire and 
Winter to arbitrary inviscid shear flow. This was then followed by the work 
of Lakshminarayana and Horlack (ref. 4) in which the effects of viscosity were 
included. These works provided the theoretical basis of our understanding of 
secondary flow phenomena. Unfortunately, except for the case of small entry 
shear or small flow deflection, these analyses do not provide a solution pro- 
cedure for general secondary flows. Researchers interested in calculating 
secondary flows have either attempted to numerically solve the governing equa- 
tions derived by these authors or numerically solve an equivalent reduced form 
of the Navier-Stokes equation (refs. 5 to 7). The present work lies midway 
between these two extremes. It combines some of the theoretical concepts sug- 
gested by the work of Hawthorne with the numerical methods suggested by the 
work of Martin (ref. 8). 

The governing differential equations which are iteratively solved in this 
paper describe a steady, compressible, inviscid and non-heat-conducting fluid 
flow in the absence of external body forces. These equations can effectively 
model many features of secondary flows (refs. 9 and 10). The compressible 
fluid is assumed to obey the perfect gas equation of state and have constant 



specific heat capacities. We will refer to the class of flows which are gov- 
erned by the above equations as Munk-Prim flows to acknowledge the applica- 
bility of the Munk-Prim substitution principle (ref. 11) which significantly 
simplifies their solution. The details of this substitution principle and its 
application to the current problem are fully discussed in Section II. Briefly, 
it entails mapping an arbitrary Munk-Prim flow into a substitute flow which is 
either homentropic (i.e., uniform entropy) or homenergetic (i.e., uniform total 
enthalpy). The choice of substitute flow is made by specifying a scalar a 
(referred to as the Munk-Prim gauge factor) which remains constant along the 
flow streamlines. A properly chosen substitute flow, as will be shown later, 
can be solved with less computational effort. Upon obtaining the substitute 
flow solution, the solution to the original problem is constructed by a mapping 
procedure. Under this mapping procedure, the pressure field, Mach number and 
streamline pattern of the substitute flow are identical to those of the origi- 
nal flow. 

The current solution procedure can also be applied to incompressible flows 
which satisfy continuity and the inviscid momentum equations, since the solu- 
tion procedure follows directly the procedures used to construct the solution 
to homoentropic compressible flows. 

The implementation of the Munk-Prim flow solver is presented on the flow 
chart shown in figure 1. The procedure is iterative and consists of two in- 
teracting loops. The derivation of the equations solved in each loop is pre- 
sented in Section I. In the outer loop, the total pressure, total temperature 
and vorticity distributions are found for a given velocity field. The solution 
strategy is presented in Section II and the questions of initial conditions 
and uniqueness of solution are addressed in Appendix A. An interesting feature 
of the current solution procedure for the vorticity vector is that it sidesteps 
solving the coupled 3-D vorticity transport equation (ref. 12). Rather the 
vorticity vector is formed from the solutions of two uncoupled partial differ- 
ential equations (PDE's). 

In the inner loop, the velocity and mass density fields. are updated so as 
to satisfy the continuity equation and the condition that the curl of the velo- 
city equals the vorticity. Within this loop, the total pressure, total temper- 
ature and vorticity distributions are assumed known. The iteration procedure 
by which the velocity and mass density fields are updated is also presented in 
Section II. This scheme is a 3-D extension of Martin's semi-direct procedure 
for solving the general 2-D Cauchy-Riemann problem (ref. 8). The boundary con- 
ditions associated with it are discussed in Appendix B while its convergence 
criterion is discussed in Appendix C. 

The results of several subsonic flow computations are presented in Section 
III to illustrate the development of secondary flow in a turning channel. For 
weak secondary flows at low Mach numbers, the computed results are shown to be 
in excellent agreement with the theoretical prediction of Squire and Winter. 

For large secondary flows two cases were analyzed, the first being the genera- 
tion of secondary flow due to inlet velocity shear, the second being the gen- 
eration of secondary flow due to inlet temperature shear. Both cases resulted 
in the generation of significant nonlinear effects which caused the flow pat- 
tern to take on a turbulent-like structure. Finally, to demonstrate the use- 
fulness of the Munk-Prim substutition principle, each Munk-Prim problem was 
solved twice, i.e., by employing a homentropic and then a homenergetic substi- 
tute flow mapping. Both procedures are shown to give identical results to 
within truncation error. 
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I - GOVERNING EQUATIONS 
( 1.1) Munk-Prim Flows 
(1.1.1) Vector Form 

The differential equations of the present analysis are the 
Continuity equation 


v • (pt) = 0 


( 1 . 1 ) 


Momentum equation 


(1/ • v)^ + — vP = 0 
p 


( 1 . 2 ) 


and Entropy equation 


1 / • vS = 0 


(1.3) 


where K, p, P, and S are the velocity vector, mass density, static pres- 
sure, and specific entropy, respectively. The above equations can be recast in 
terms of total pressure P Q , total temperature T q , and velocity 1 / by employ- 
ing the following thermodynamic relationships: 


p 




1 



2 C T I 
P o) 


(1.4) 


P = 


■\i . n 

2 rrr 

P O) 



and 


S = C In T - R In P 
p o o 


(1.5) 


( 1 . 6 ) 


Here C p is the specific heat capacity at constant pressure, R the ideal gas 
constant, and y the ratio of specific heat capacities. The parameters C , 

r 

R, and y are all constants for Munk-Prim flows. Throughout this study, the 

focus is to find the five independent variables P , T Q , and 1 'I which satisfy 

the five independent equations (eqs. (1.1) to (1.3)). All other variables are 
considered functions of these fundamental variables. 
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To simplify the solution of the numerical problem, equations (1.1) to (1-6) 
are nondimensional ized by scaling length, pressure, temperature and density 
against a reference length L*, a reference pressure P*, a reference tempera- 
ture T*, and a reference density p*(P*, T*, and p* are assumed to satisfy 
the perfect gas law). Any additional reference variable is formed directly 

from P*, T*, and p*. (For example, the dimensionless velocity d =^* - 

and the dimensionless entropy d = f * — - ^ with S* d = f * C, 

*' I 

This results in a set of dimensionless flow equations which retain their ori- 
ginal form with a dimensionless gas constant R = 1 and a dimensionless heat 

capacity C p = ^ 

In the remainder of this subsection, we will collect the differential 
equations which are used in the two iteration loops which form the solver. To 
proceed, one notes that equation (1.2) and the thermodynamic relations can be 
combined to form Crocco’s equation, i.e. 


y/WT^* 

In T* - R In P*.) 


1 x ft = vh Q - T vS 


(1.7) 


Here h Q and T are the specific total enthalpy and absolute flow temper- 
ature, respectively. The vorticity vector a is defined as 

Q = v x 17 (1.8) 


and since v • (v x H) = 0 where is an arbitrary differentiable vector 
field, the vorticity field must be solenoidal, i.e., 

v • n = 0 (1.9) 


In the current analysis Crocco's equation is used in place of the momentum 
equation. Also, since v • ! x ft = 0, equations (1.3) and (1.7) imply that 

17 • v h Q = 0 (I. 10) 


i.e., the specific total enthalpy is conserved along a streamline. With the 
aid of equation (1.6) and the relation: 


h 


o 



(I. ID 


Equations (1.3) and (I. 10) can be replaced with the equivalent expressions: 


17 • vP Q = 0 


( 1 . 12 ) 


and 
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(1.13) 


t • vT 0 = 0 


In the inner iteration loop, P , T , and n are provided by the outer 

loop and the objective is to find a density p and a velocity 1 1 which 
satisfy equations (1.1) and (1.8). The details of the inner loop iteration 
procedure are provided in Section II. 

In the outer loop, the velocity field D is assumed known and the objec- 
tive is to find a total pressure P Q , a total temperature T , and a vorticity 

ft which satisfy equations (1.7), (1.9), (1.12), and (1.13). The evaluation of 
these variables can be divided into three steps. First, the unknowns P Q and 

T q are solved using equations (1.12) and (1.13). This can be done if the in- 
let profiles of P and T are known. In the second step, the variables 
oo 

h Q , T and S which appear on the right side of equation (1.7) are evaluated 

from the updated P Q and T q field. In the final step, the vorticity vector 

ft, which is required to satisfy equations (1.7) and (1.9), is constructed 
according to Theorem 2 of Appendix A. The .result is 


-► 

fl = 

where the two scalar functions 
and 


vh Q x Vt + vS x vu 

(1.14) 

t and g satisfy the equations 


1 1 • Vt = 1 

(1.15) 

• vu = -T 

(1-16) 


(In deriving eqs. (1.14) to (1-16), it was assumed that "V ^ 0 and either 
vh Q £ 0 or vS ^ 0 within the solution domain.) The motivation for express- 
ing ft in the form of equation (1.14) is that it can be obtained directly from 
the solution of two uncoupled first order PDE's. This procedure is far sim- 
pler than trying to solve equations (1.7) and (1.9) directly. The solution of 
equations (1.15) and (1.16) requires as input the inlet profile of t and p. 
These profiles, according to Theorem 1 of Appendix A, must be chosen to be con- 
sistent with the streamwise vorticity at the inlet. 

In the present analysis, the substitute flows are either homentropic 
(vS = 0) or homenergetic (vh Q = 0). For homentropic flows, the second term on 

the right side of equation (1.14) vanishes. As a result, there is no need to 
solve equation (1.16). Furthermore, for this class of flows, T 0 becomes a 
function of P 0 (see eq. (II. 5)) which eliminates the need to solve equation 
(1.13). Thus only two PDE's, i.e., equations (1.12) and (1.15) need to be 
solved in the outer loop for homentropic flows. Similarly, for homenergetic 
flows, only equations (1.12) and (1.16) need to be solved. This observation 
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justifies our earlier contention that the use of the Munk-Prim substitution 
principle can lead to a relatively simple procedure for solving the general 
Munk-Prim problem. 


(1.1.2) Tensor Form 


In the previous subsection, the basic equations for the two iteration 
loops were derived using vector notation. Before these equations are solved in 
a non-Cartesian computational space, they shall be converted to general tensor 
forms. If we adopt Einstein's summation convention (which we shall use 
throughout this paper), then in terms of the general computational coordinates 
{x 1 }, the tensor forms of the outer loop equations (1.12) to (1.16) are, re- 
spectively. 


and 


, aP. 
f 1 — 2- 
ax 1 


ax’ 


= 0 


= 0 


( 3 h 

o^ + aS ^ 

ax” ax ax J ax 


V 1 1 

ax 1 


V 1 3y — _ _t 

ax 1 


(1.17) 

(1.18) 

(1.19) 

( 1 . 20 ) 

( 1 . 21 ) 


In the above equations V 1 is the contravariant velocity vector, a 1 is 
the contravariant vorticity vector, g is the determinant of the covariant 

i ik . 

metric tensor g . . , and e J is the contravariant Levi-Civita tensor density. 

* J 

It should be understood that the range of tensor indices, unless otherwise 
noted, is 1, 2, 3. 

The tensor forms of the inner loop equations (eqs. ( 1 . 1 ) and (1.8)) are, 
respectively, 

aF^ 

= 0 ( 1 . 22 ) 

< ax' 


and 
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(1.23) 


\/g fl 1 = e 


ijk 


3V, 


3X 


where F 1 is the contravariant mass flux vector density defined as 


-i def. 


pV? 9 1J Vj 


(1.24) 


The variables which appear in equations (1.23) and (1.24) and were not defined 
previously are V-j, the covariant velocity vector and g 1 J the contra- 
variant metric tensor. 


(1.2) Incompressible Flows 

The analysis of the previous subsection can be extended to incompressible 
flows. For these flows, the governing equations are equations (1.1), (1.2), 
and 


p = constant (1.25) 

These equations can be nondimensional ized by introducing a reference length 

d 

L*, a reference density p*( = * fluid density), a reference velocity V* and 

a reference pressure P*( d = f * p * (V*) 2 ). This results in a set of dimension- 
less flow equations which retain their original form with a dimensionless den- 
sity of p = 1. 

Combining equations (1.2), (1.8), and (1.25) yields 

txsi = v(P o /p) (1.26) 


with 


p ^ef. p + 
o 



= total pressure 


(1.27) 


As a result of equations (1.25) and (1.26), equation (1.12) is also valid for 
this special class of flows. 

Equations (1.1), (1.8), (1.9), (1.12), and (1.26) are the equations to be 
solved for incompressible flows. Their solution can be obtained in a fashion 
similar to that for the homentropic flows. 


II - SOLUTION STRATEGY 

1 1 . 1 Munk-Prim Substitution Principle and Its Application 

In Section I, it was shown that a Munk-Prim flow can be specified in terms 
of total pressure P , total temperature T q and the velocity vector D. The 


7 



Munk-Prim substitution principle states that, if P , T , and t specify a 

given Munk-Prim flow, then one can generate another Munk-Prim flow (referred to 
as the substitute flow) by the following transformation: 

p 0 - p o ! T 0 - V“ 2 ' v =1// “ (n.1) 


where PJ, and "V 1 are the flow variables of substitute flow and o (Munk- 
Prim gauge factor) satisfies the equation: 

• 7 a = V • va = 0 (II. 2) 

The proof of this principle follows directly from equations (1.1) to (1.6). 

From equation (1.5) and equation (II. 1), it is seen that the static pressure is 
invariant under a Munk-Prim transformation. Similarly, the Mach number and the 
streamline pattern are also invariant under this transformation. 

The substitution principle associated with Munk-Prim flows can be used to 
reduce the computational effort required to solve the general Munk-Prim flow 
problem. This reduction in effort is made possible by a proper choice of a. 

To explore this further, let a be equal to 


a 




(II. 3) 


or 


a 




(II. 4) 


Both a. and ct 0 satisfy equation (II. 2) since T and P are constants 
i i o o 

along a streamline. If a = then equations (II. 1) and (II. 3) coupled with 
equation (1.6) imply that 

Yzi xd 

S' = 0 and T' = P' y = P Y (II. 5) 

0 0 0 


i.e., the substitution flow is homentropic. Similarly, if a = c^, 

T* = 1 and P' = P n 
0 0 0 - 


(II. 6) 


i.e., the substitute flow is homenergetic. Solutions with either S' = 0 or 
Tq = 1 require less work to evaluate since their generation requires two 
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fewer equations to be solved per outer iteration cycle than that for the gen- 
eral flow problem. 

The inlet conditions of a substitute flow can be evaluated from those of 
the original flow with the aid of equations (II. 1), (II. 3), and (II. 4). Using 
these conditions and the solid wall boundary conditions (these conditions are 
invariant under the Munk-Prim transformation) along with a specified outflow 
condition, a substitute flow solution can be computed. 

To generate the solution to the original problem from the solution of the 
substitute flow problem, the Munk-Prim gauge factor ^(c^) must be evaluated 

throughout the flow field. This is accomplished by integrating equation 
(II. 2). With c^, ( 012 ) known, the solution to the original problem is found by 

employing equation (II. 1). 


1 1.2 Outer Loop 

In Section I, it was shown that only equations (1.17) and (1.20) or equa- 
tions (1.17) and (1.21) need to be solved in the outer loop if a solution to a 
substitute flow problem is desired. These PDE's are first order and linear if 
the velocity field is assumed known. They may be solved using a standard 
marching procedure for hyperbolic equations. .The algorithm used to integrate 
these equations in the present study was constructed based on the work of 
Hartree (ref. 13). It will be explained using equation (1.20) as an example. 

We assume that the computational domain is a paral lelpiped which is bounded by 

2 3 

four solid walls in x and x directions (see fig. 2(a)). We further as- 
sume that V 1 is positive everywhere and that all the values of t are known 
on the grid plane GP1. 

For an interior grid point A on GP2 at which t is to be evaluated, a 
point Q' on GP1 is found (see fig. 2(b)) such that the vector Q'A is in the 
direction of the known velocity vector at A (denoted by V n (A)). The velocity 

vector at Q' (denoted by V 1 ( Q * ) ) is computed by linearly interpolating (four 
point bivariate interpolation (ref. 14)) the velocity vectors at the four sur- 
rounding grid points A', B', C' and O'. The average of V 1 ( Q 1 ) and V n ( A) is 
then evaluated and set equal to V n (A). Next a point Q is located on GP1 such 

that the vector QA points in the direction of V^A). The value of T at 
point Q (denoted by x(Q)) is then found by linear interpolation using the 
known value of t at A', B 1 , C', and D'. In the final step, t is evaluated 
at point A (denoted by t(A)) from the equation 

t(A) = t (Q) (II. 7) 

V \A) 


where <sx is the distance between GP1 and GP2. If point A is a boundary grid 
point, the above procedure becomes two-dimensional (the normal velocity compo- 
nent vanishes at the channel walls). The integration of the three remaining 
transport equations (eqs. (1.17), (1.18), and (1.21)) can be performed in a 
similar fashion. 
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Once the transport equations are solved, the vorticity field is updated 
according to equation (1.19). The finite difference form of this equation is 
obtained by approximating grad with a second order central difference operator. 
The only exception occurs at. the boundary where the normal derivative which 
appears in grad is approximated by either a forward or backward second order 
difference operator. 


1 1.3 Inner Iteration Loop 

The basic equations which are solved in the inner loop were identified in 
the previous section. They are equations (1.22), (1.23), and (1.24) and are 
repeated here for clarity. 


a F 1 

= 0 (II. 8) 

3 X ' 


aV 

\/g a 1 = e ijk — j (II. 9) 

3X J 


F 1 = P yg g ij V ‘ (II. 10) 

J 

where, for Munk-Prim flows, the density is related to the velocity by equation 
(1.4). This equation, expressed in general tensor form, is 

1 


p 



s’VjV 

2 s v 


(II. 11) 


One is reminded that P Q , T , and n 1 are known variables in the inner loop. 

To solve equations (II. 8) to ( I I— 10) , we introduce the following iterative 
scheme which is suggested by the work of Martin (ref. 8): (n = 1, 2, 3,...) 


SF 1 '") ^ 8 F 2 <"> „ 8 F~ 3(n) 
1 — 2 — 3 

3 X 3X 3X 


( 11 . 12 ) 


qP 3 ( n ) a p2(n) aF 3(n-l) aF 2(n-l) 

ax 2 ax 3 3X 2 3x 3 

aF^ n ^ 3F 3(n) _ aF^ n ~^ aF 3 ^ 11 " 1 ^ 

3X 3 3X 3X 3 3X^ 


(11.13) 


(11.14) 
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(11.15) 


p2(n) 3 pl(n) a p2(n-l) 3 pl(n-l) 


3X 1 

3X 2 

3X 1 

ax 2 

avt"> 

3 vl n ^ 

+ "5 

3vJ n_1) 

aV^-D 

■4- £ 

3X 2 

ax 3 

3X 1 

3X 2 


,v("> 

3V^ n) 

H 

C! 

*> 


ax 2 

ax 3 ' 


aV<"' 

>n n) 

OJ 

cs 

§ 


ax 3 

ax 1 " 


aV< n) 

av‘"> 



ax 1 

1 

C\J 

X 

cx> 


(11.16) 


(11.17) 


(11.18) 


(11.19) 


Here F 1 and are, respectively, the values of F 1 and V 1| at 

the end of (n-l)-th iteration. F^ n ^ and are, respectively, the inter- 

mediate values of F 1 and V at n-th iteration. The above iteration scheme 

is motivated oy reasons which will be made clear shortly. For now, it suffices 
to note that these equations uniquely define all of the intermediate variables 
at n-th iteration if all of the variables at the end of (n-l)-th iteration are 
known and proper boundary conditions are imposed. This iterative scheme is used 
solely in a given computation space, therefore, equations (11.12) to (11-19) and 
any equations to be derived from them shall not be considered as tensor equa- 
tions. However, in order to economize the presentation of the iterative equa- 
tions to follow, we shall continue to use tensor notations and Einstein conven- 
tion with one important exception, i.e., the usual tensor analysis convention 
concerning the positioning (upper or lower) of indices will be disregarded . 

Equations (II. 12 ) to (11-19) do not provide the mechanism tor updating the 
variables without tilda. Six additional equations must be supplied to perform 
this task. Two of the many possible sets of relationship are given by the 
equations 

Scheme 1: = vj n ^ (11.20) 

F lln) = P (n) yg g 1j V< n > ' (11.21) 
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Scheme 2: 


F i (n) = pi (n) 


( 11 . 22 ) 



7 ^ 





(II 23) 


where is related to by equation (II. 11) for Munk-Prim flow. 

If Scheme 1 is used to update the variables without tilda, then and 

become a function of vj n ^. The tilda variable F^ n ^ plays no role in 


the updating of vj n ^ and F| n ^ and hence only equations (11.16) to (11.19) 

need to be solved. Conversely, if Scheme 2 is used, then both F 1 ^ and 

vj n) become functions of F^ n ^ and only equations (11.12) to (11.15) need to 

be solved. It is noted that the use of Scheme 1, effectively solves equation 
(II. 9) while leaving equation (II. 8 ) unsolved. This fact follows directly from 
equations (11.17) to (11.20). Conversely, the use of Scheme 2 solves equation 
(II. 8 ) while leaving equation (II. 9) unsolved. For the current solver. Schemes 
1 and 2 are used alternatively. In other words. Scheme 1 is used initially 
and at all subsequent iterations for which n is odd while Scheme 2 is used 
when n is even. If this iterative procedure converges, the resulting solu- 
tion must satisfy equations (II. 8 ) to (II. 11). 

To simplify equations (11.12) to (11.19), we first introduce the variables 

f^ n ^ which are defined by: 


f i (n) def. pi(n) _ pi(n-l) 


(11.24) 


Equations (11.13) to (11.15) coupled with equation (11.24) imply that 


ijk 3f 


k(n) 


= 0 




(11.25) 


As a result, there exists a scalar function 9 ^ such that 

f i(n) _ a<p ^ 

3X 1 

Combining equations (11.24) and (11.26), one obtains 

pi (n) = pi(n-l) + 

ax 1 


(11.26) 


(11.27) 


Substituting equation (11.27) into equation (11.12), the Poisson's equation for 
i.e., 


12 



(11.28) 


a 2 ,("> 3 F i(n - 1} 

ax 1 ax 1 ax 1 


is obtained. Thus, we have successfully reduced equations (11.12) to (11.15) 

to a single Poisson's equation (eq. (11.28)). If q>^ is known, F 1 ^ can 
be calculated using equation (11.27). 

Next equations (11.16) to (11.19) are transformed into a set of Poisson 

n ) 

equations. We again introduce a correction variable f^ which this time is 
defined as: 

x (n) def. y(n) _ v (n-l) 
i i i 

Using equation (II 29), equation (11.16) becomes 

' r\ 


(11.29) 


(11.30) 


Thus there is a vector potential 



such that (ref. 15) 





3X J 


(11.31) 


and 




(n) 


= 0 


3 X 


Combining equations (11.29) and (11.31), one obtains 


w( n) v (n-l) + ijk 341 k 
i i 

3X V 


(n) 


(11.32) 


(11.33) 


Substituting equation (11.33) into equations (11.17) to (11.19) and using 
equation (11.32), we obtain 


.2,00 

3 41^ 

3X J 3X J 


= e 


3 V ( n_1 ) 

ijk 3V k 

3X^ 


- V9 


(11.34) 
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Equations (11.16) to (11.19) have thus been reduced to three Poisson's equa- 

v(") 

(n) i 

in terms of ^ . At first glance, equations (11.32) and (11.34) give the 


tions (11.34) and a guage condition (11.32). Equation (11.33) defines 

equa 
(n). 


impression that the three unknowns ^ 's are overspecified by four indepen- 
dent differential equations. This, however, is incorrect as will be shown. 

We first noted that the tensor form of equation (1.9) is 


a(vW) 


= 0 


3 X 


(11.35) 


which coupled with equation (11.34) yields 


3X J 3X J 



= 0 


(11.36) 


Thus 



satisfies Laplace's equation if satisfies equation (11.34). 


If we choose the boundary conditions (B.C.S) for ^ 


(n) 


such that 



vanishes at all the boundaries of the computational domain, equation (11.32) 
is the unique solution of equation (11.36). In other words, the gauge condi- 
tion (11.32) is satisfied in the entire solution domain if and only if it is 
satisfied at the boundaries. 

Equations (11.16) to (11.19) can be simplified in yet another fashion if a 


k of equation (II. 9) is known. 


particular solution U 
definition, is required to satisfy the equation 


The vector U k , by 


A - V ^ 1 

3X J 


(11.37) 


If, in addition, we define a variable f^ n ^ such that 

f-(n) = y(n) _ u _ (11.38) 

Then equations (11.17) to (11.19) and (11.37) imply that the Curl of f^ 
vanishes, i.e.. 


3f^ n) 

e ijk -X- = 0 
3X J 


(11.39) 
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Then, there must exist a scalar such that 

7 (n) _ a q (n) 

1 ax 1 


(11.40) 


which when combined with equation (11.38), yields 

(n) 


V< n > = II + 12 — 

1 ' 1 ax 1 


(11.41) 


Upon substituting this result into equation (11.16) one obtains 


a 2 ( n ) 

3 a 

ax 1 ax 1 


3 

ax 


T fS"- 1 ’ - ",) 


(11.42) 


This time, equations (11.16) to (11.19) have been reduced to a single Poisson's 

equation (eqs. (11.42) and (11.41)) which updates in terms of U. 

and 


' Of the two methods presented for solving equations (11.16) to (11.19), the 
second method has the advantage of solving only one Poisson's equation (eq. 
(11.42)) instead of three (eq. (11.34)). Its disadvantage is that it requires 
a particular solution of U-j of equation (II. 9) to be known. To circumvent 
this problem, the following strategy is adopted: Each time through the inner 

loop, the vector potential method (i.e., the first method discussed) is used to 

obtain . Assuming = vj^, the second method is then used to obtain 

y(3), v^,... at subsequent odd iteration cycles. Since and F^ n ^ 


are updated in terms of if n is odd and in terms of F 1 ^ if n 

1 ~i(2) 

is even, the above strategy combined with the method for updating F v , 

~ i ( 4 ) 

F ' .... leads to a complete iterative procedure for solving equations (II. 8) 

and ( 1 1 . 9) . This procedure, except for the first iteration cycle, involves 
only recursively solving two Poisson's equations, i.e., equation (11.28) for 
even n cycles and equation (11.42) for odd n cycles. The construction of the 
iteration procedure equations (11.12) to (11.19) made this possible. 

The boundary conditions required to solve the Poisson's equations are 
presented in Appendix B. The stability of this scheme is discussed in Appen- 
dix C and shown to be a function of the metric tensor g^.. 

The equations which appear in the inner loop will be solved numerically. 
The divergence, grad and curl are approximated by second order central differ- 
ence operators in the interior of the computational domain. At the boundaries 
the normal derivative operator which appears in these expressions is approxi- 
mated by either a second order forward or backward finite difference operator. 
The finite difference form of the Laplacian is fixed by the Fast Solver (ref. 
16) which we use to solve the Poisson equations. 
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Ill - APPLICATIONS TO A SHEAR FLOW IN A TURNING CHANNEL 


III.l Problem Definition and Solution Procedure 

In this section, the solution procedure outlined in the previous section 
is used to compute the developing inviscid rotational flow field in a turning 
channel (fig. 3). This channel is generated from a parallelpiped (fig. 4) in 
computational space using the following transformation: 


e uX = J ^COSh(nX^) - cos(nx 2 )^ 

(III.l) 

irx 2 1/2* Sinh(irX 1 /2)C0S(irX 2 /2) 

COS q — 

(III. 2) 

|/c0Sh(irX 1 ) - C0s(nx 2 ) 

-3 3 

X = X 

(III. 3) 


—1—2—3 12 3 

where the coordinates (x , x , x ) refer to physical space and (x , x , x ) to 

i o i o 

computational space. The transformation between (x , x ) and (x , x ) is a 
conformal transformation defined by 

Z = (2/ir)ln(sinh(irW/2)) (III. 4) 


— 1—2 12 

with Z = (x , x ) and W = (x ,x ). The parameters a, b, c, and d which de- 
fine the flow boundaries are subjected to the conditions -» < a < b < +® and 

3 

and 0 < c < d < 1. The floor and ceiling of a channel are at x = 0 and 

3 

x = e, (0 < e < -«*>). According to equatuions (III.l) to (III. 3), the values 
of the metric tensor in computational space is given by 



(III. 5) 


where 


1 2 

cosh(Tix ) + cos( ttx ) 

1 2 

C0Sh( irX ) - C0S(irX ) 


(III. 6) 


Equation (III. 5) implies that two orthogonal coordinate lines in computational 
space remain orthogonal in physical space. Also according to equation (III. 5) 

and (III. 6), g • • -*• 6.. (Kronecker Delta symbol) as |x^| + +®. Therefore, in a 

* J ‘ J 
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region where x* is sufficiently large, the flow equations, the flow descrip- 
tion and the channel geometry as viewed from both computational and physical 
spaces are identical. We shall assume that both the inlet and exit planes lie 
in such a region. 

The first step in solving a Munk-Prim flow problem, as seen from figure 1, 
is the specification of the inlet flow conditions. The allowable inlet condi- 
tions will be limited to those of a parallel entrance flow. Hence, according 
to equations (8.1) to (B.3), the inlet conditions are specified by an arbitrary 
axial velocity profile, an arbitrary temperature profile and a uniform pressure 
profile. In the current study, the inlet conditions are explicitly given by: 
(Since v ~ 1 at the inlet, these conditions can be given directly in computa- 
tional space) 


v 1 = v c • u(x 3 , 6 V ), V 2 = V 3 = 0 


(III. 7) 


T = 




(III. 8) 


P = 


(III. 9) 


Here V , T , P , 6 , and a T are constants and 

^ L « v I 



a) is defined by 


u(x 3 , a) = 1 + a • 



(III. 10) 


where e is the height of the channel. For a fixed value of the constant a, 

3 

u is a normalized linear function of x with its maximum variation given by 
a , i . e . , 


1 

e 




(III. 11) 


and 


a = u(e,a) - u(o,a) 


(III. 12) 


As a result of this normalization, the constants V c and T c are the average 
inlet axial velocity and the average inlet temperature. The constants P and 

V 

T c always can be assigned a value of 

P c = T c = 1 (III. 13) 
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by a proper choice of the reference state variables P* and T*. The inlet 
flow conditions can thus be specified solely in terms of V , 5 , and 6 T . , 

v# V I 

From these conditions, the inlet flow profiles of a substitute flow can be de- 
rived using equations (B.4), (B.5), (II. 1), (II. 3), and (II. 4). The inlet pro- 
files of the function t and y (see eq. ((1.14)) which are consistent with 
the above inlet flow conditions are given by equation (B.6). 

Prior to the start of iteration the substitute flow is initialized so that 

the density p and the contravariant mass flux vector density F 1 satisfy the 
following conditions: 


= 0, = 0, F 2 = F 3 = 0 

ax 1 ax 1 


(III. 14) 


These equations coupled with the given inlet profiles for p and F* uniquely 
determine the initial values of p and F 1 within the computational domain. 
Subsequently, the contravariant velocity field V 1 can be evaluated through 

the relation: V 1 = F 1 /(pVo)‘ The flow field so initialized is consistent 

with the inlet flow conditions, the solid wall boundary conditions (eq. (B.8)) 
and the continuity equation (eq. (II.8) ) . 

Following the flow initiation, the solver will iterate between the inner 
and outer loops until the values of P q , T q , and n generated by two consecu- 
tive outer loops converge to a preassigned tolerance level. The converged so- 
lution to the substitute flow problem is then used to construct the solution 
to the original problem according to the theory of Section II. 1. 

The finite difference solutions to the various equations which form the 
inner and outer loops are constructed using the basic grid shown in figure 5. 
This grid is channel conforming with uniform spacing in each direction in com- 
putational space. This uniformity, however, is lost in physical space. 

The stability of the algorithm for solving the transport equations is in- 
sured by subdividing each basic grid interval ax^ into M equally spaced 

subintervals. The length of these subintervals, 6x\ is chosen such that the 
stability criteria 


sx 


< 1 and 


6X ' 


ax 


< 1 


(III. 15) 


2 3 

is satisfied everywhere on the refined grid. In this equation ax and ax 

2 3 

are, respecti vely, the lengths of the basic grid intervals in the x and x 
directions. The refined grid velocity field is obtained by linearly interpo- 
lating the basic grid velocity field provided by the inner loop. 
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III. 2 Numerical Results i 

' * * * f . 

This subsection will be devoted to assessing the ability of the solver to 
capture the physics of the development of inviscid secondary flows in a turning 
channel. In order to perform this assessment, six problems were analyzed. For 
each problem, the values of the parameters previously defined are tabulated in 
Table I. In addition, the values of the grid parameters NX1, NX2, and NX3 

12 3 

which define the number of grid intervals in the x , x , and x directions 
are also tabulated along with the value of NX1S which equals the number of sub- 

intervals in a grid interval ax^. The value of the parameters c and d were 
chosen to yield a channel of nearly constant width (approximately equal to 
(d-c)) in physical space. For all six problems the specified inlet flow condi- 
tions and channel geometry insure subsonic flow conditions. The inlet Mach 

number based on V and T (=1) for the first two cases is 0.0845: for the 

c c 

remaining four problems, its value is 0.4. 

The first two problems were constructed to examine the ability of the sol- 
ver to reproduce the theoretical results of incompressible weak secondary flow 
theory. For these problems the specified inlet flow conditions lead to the 
establishment of a homentropic flow. Thus these flows were calculated without 
the aid of the Munk-Prim substitution principle. For the limiting case of in- 
compressible weak secondary flow, Squire and Winter (ref. 1) derived the 
expression 


5 - 2 v 


(III. 16) 


for the growth of secondary vorticity along a streamline. In this equation, n Q 

is the inlet vorticity of the streamline, £ is the secondary vorticity at any 
point on the streamline where the local angle of deflection (i.e., the angle 
between the streamline direction at the inlet and its direction at a point in 
the channel] is e. Equation (III. 16) implies that, along a streamline, the 
secondary vorticity increases linearly with 0 . 

For the chosen channel the values of e are more or less constant across 
any cross section. If the imposed inlet flow field is such that all stream- 
lines share a common value of ft (n = 6 /e for the current problems) it is 

o' 0 v 

seen that £ is nearly uniform across any cross section of the channel. The 
results to be presented for stations A, B, C, D, and E of figures 1 and 2 will 
be shown to be consistent with these observations. These stations are channel 

cross sections defined by x^ = -1/3, x* =0, x^ = 1/3, x* = 2/3, and x* = 1, 
respectively. In physical space these. surfaces are slightly curved. In pre- 
senting our results in physical space, we, neglect the effect of surface curva- 
ture and consider these stations to be flat. .This approximation does not in- 
troduce any significant errors. 

• The secondary velocity fields at Stations B and C for the first problem 
are shown in figure 6. Comparing the magnitude of the displayed velocity vec- 
tors to the inlet area avaraged velocity of V = 0.1 (the -length scale of the 

vectors is indicated by a standard vector of magnitude 0.001 at the bottom- 
right corner of each plot), they are seen to be very small. Thus the theory 
of Squire and Winter should be applicable to this problem. Both figures show 
the generated secondary velocity field exhibits a sol id' body rotational pattern 


19 


over most of the cross section, which indicates the secondary vorticity is 
nearly uniform over each cross section. The only exception occurs in the cor- 
ner regions where equation {III. 16) ceases to be applicable due to excessive 
streamline curvature. 

Furthermore, figure 6 indicates that one may view the line defined by 

2 3 

x = 0.5 and x = 0.05 (according to fig. 4, this line passes through the 
central grid point at all cross sections) as a streamline since the secondary 
velocity is nearly zero along this line. At the five points where this approx- 
imate streamline intersects Stations A, B, C, D, and E, the values of e, £, 
and u/(2 n Q e) for the first two problems are tabulated in Table II. The nu- 
merical results are seen to be in excellent agreement with equation (III. 16). 
These results also indicate that, as the flow proceeds further down the chan- 
nel and the secondary velocity grows larger, nonlinear effects begin to come 
into play and the agreement between linear theory and numerical results 
deteriorates. 

For the third problem, a significantly larger value of inlet velocity 
shear was chosen to assess the ability of the solver to capture the physics 
associated with large secondary flows. Figure 7 shows the secondary velocities 
that developed at Stations B, C, and E. At Station B, the secondary velocity 
field is well established, exhibiting the solid body rotation pattern predicted 
by the linear analysis of Squire and Winter. As the flow proceeds further 
around the bend, the local secondary velocity field near the four corners grows 
in magnitude, with a maximum approaching the average inlet velocity V c = 0.4733 

(a vector whose magnitude equals the average inlet velocity of 0.4 is depicted 
at the lower right corner of each figure). This flow structure can no longer 
be predicted by the linear theory of Squire and Winter. At station E which 
lies well downstream of the channel bend, the nonlinear behavior grows more 
evident as two additional vortices appear in the upper-left and lower-right 
corners. Further down the channel, the secondary flow field becomes more ill 
defined and begins to take on the appearance of turbulence. 

To study other aspects of this homentropic flow, the contours of constant 
total pressure P Q at station B, C, and E were computed and are shown in fig- 
ure 8. Since a surface of constant P Q is a stream surface (see eq. (1.12)), 

these illustrations also represent the cross-sectional views of several stream 
surfaces. For the current problem, the inlet flow conditions vary only in the 

3 

x -direction. Therefore, the stream surfaces are horizontal at the channel 
entrance. The contour values of P Q were chosen to yield contours which coin- 
cide with the horizontal grid lines at the entrance to the channel. 

As the flow proceeds down the channel, the stream surfaces are rotated by 
the secondary velocity field to near vertical slightly downstream of Station B. 
The only exception occurs at the upper-right and lower-left corner regions 
where they tend to become clustered. This phenomenon can be explained by re- 
calling that every stream surface depicted in figure 8 is swept out by a group 
of fluid particles which originally passed through a horizontal grid line at 
the channel entrance. As the flow proceeds down the chartnel, the fluid parti- 
cles initially attached to the right (left) side wall will remain attached as 
long as streamline bifurcation does not take place. In addition these parti- 
cles are transported toward the top or bottom walls by the secondary velocity 
field. Since the top and bottom walls are stream surfaces on which P is 

either a maximum or a minimum (i.e., maximum on the top wall, minimum on the 
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bottom wall), the fluid particles on the side walls can never reach the top or 
bottom unless equation (1.12) is violated. They become stagnated when viewed 
in the cross-sectional plane. This leads to the clustering phenomena depicted 
in the figure. At Station E, the total pressure contour becomes highly con- 
torted and also shows additional clustering in the upper-left and lower-right 
corner regions. 

The physics of the flow field illustrated in figure 8 is very complex and 
highly nonlinear. To estimate the accuracy of this solution and, particularly, 
the impact on solution accuracy of the first order integration scheme for the 
transport equation, a grid refinement study was undertaken. This study con- 
sists of Problems (3), (4), and (5). These three problems represent the solu- 
tion of the same flow problem (all have identical values of geometric and flow 
parameters) on three different grids whose cells (in computational space) are 
geometrically similar. The cell length ratio as well as the subinterval length 
ratio is 6:4:3. This construction insures that any numerical difference be- 
tween the solutions to these problems is characterized by a single parameter, 
i.e., the cell length in any one direction. 

Figure 9 shows the contours of constant total pressure computed on the 
coarse grid at Stations B, C, and E. The contour values of total pressure are 
again chosen to yield contours which coincide with the horizontal grid lines 
at the channel entrance. Superimposed on these plots are the corresponding 
total pressure contours computed using the intermediate and finest grids. At 
station B (half way around the channel), the results for all three grids are 
in close agreement. The slight differences one observes are confined to the 
channel side walls and can be attributed to the inadequate resolution of the 
clustering of stream surfaces in the corner regions. 

As the flow proceeds downstream, the influence of these clustered stream 
surfaces on the accuracy of the results grows as evident by the contour plots 
for Station C. These plots show that the results for the coarse and fine grids 
are only in qualitataive agreement, while the results for the intermediate and 
fine grids still remain in close agreement in the central portion of the 
channel . 

On the side walls, nearly all the contours have been compressed into a 
single grid interval. The accuracy of the results in this region are therefore 
questionable and will not be addressed any further in this study. The final 
comparison is for Station E. The agreement between the three sets of results 
has suffered additional deterioration, with the results for the intermediate 
and finest grids now showing only qualitative agreement. One concludes from 
this study that, except in the corner regions, the results presented in fig- 
ures 8 are independent of grid size for Stations B and C and can be further 
improved with additional grid refinement for Station E. 

The secondary flows which we have discussed so far have all been generated 
as a result of an inlet velocity shear. In Problem (6), we examine the second- 
ary flow as generated by a temperature shear. The entering velocity field is 

3 

uniform, while the inlet temperature and entropy vary in only the x direc- 
tion. The solution to this general Munk-Prim problem was obtained assuming the 
substitute flow to be homentropic. 

The resulting secondary velocity fields developed at Stations B, C, and E 
are shown in figure 10. Comparing these results with those on figure 7, it is 
seen that the secondary velocity fields are very similar except for their sense 
of rotation. This unexpected result can be explained by recalling that the 
streamline pattern for Problem (6) is identical to that of the homentropic sub- 
stitute flow from which it was constructed. This substitute flow has a nega- 
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3 

tive velocity shear (i.e., axial velocity decreases in the x direction) and 
a uniform temperature and entropy profile at the channel inlet. These inlet 
flow conditions with the exception of the sign of the inlet velocity shear are 
very similar to those for Problem (3). This accounts for their close resem- 
blance. Furthermore, examination of the contours of constant total pressure 
(not shown) further corroborated the close resemblance of these two problems. 

The flow problems discussed thus far have either been solved without the 
direct use of the Munk-Prim substitution principle (Problems (1) to (5)) or 
solved assuming a homentropic substitute flow (Problem (6)). They could have 
been just as easily solved assuming a homenergetic substitute flow. Numerical 
results based on this alternative formulation for all six problems verified 
this contention. In fact, the graphs (either the secondary velocity fields or 
the constant total contours) generated by these two approaches were indistin- 
guishable. This result, along with the successful computation of Problem (6), 
demonstrates the reliability of the Munk-Prim substitution principle as a tool 
for steady inviscid flow computation. 


CONCLUSIONS 

The current algorithm represents a new procedure for solving inviscid sub- 
sonic 3-D rotational flows. The convergence for the algorithm is shown to be 
dependent on the value of the matrix tensor in computational space. Results 
obtained using this algorithm show that it can accurately predict the behavior 
of weak secondary flows. Furthermore, the results of a grid refinement study 
indicate the algorithm is sufficiently accurate to capture the physics of the 
evolving strong shear flows. 

Finally, this work clearly shows that the use of the Munk-Prim substitu- 
tion principle can substantially reduce the computational effort n solving 
any 3-D inviscid steady rotational flow problem. 
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APPENDIX A 


Basic Theorems of the Outer Loop 

In the current study, the outer loop solves for the vorticity field n 
which satisfies the equations 


1 / x a = vh Q - T vS 


(A. 1) 


▼ A 
v • a = 0 


(A. 2) 


Here, "V, h , T, and S are assumed to be known flow variables with h Q and S 
satisfying the equations 


1/ • vh Q = 0 


(A. 3) 


"V • "v S = 0 


(A.4) 


Given the above conditions, we establish the following theorems: 


Theorem 1: Assuming 1/^0 throughout the region of interest and defining 

~\j 

e v = -j — r then a solution n of equations (A.l) and (A. 2) is 


uniquely defined if the streamwise vorticity (n • e ) is specified 
at the inlet surface. 

Proof: Let n' and a 11 be two solutions of equations (A.l) and (A. 2) 

with n 1 • e v = n" • e at the inlet surface. If we define 

= n' 1 - n ' , then one has 


ll x 6n = 0 


(A. 5) 


v • = 0 


(A. 6) 


and, at the inlet surface. 


<sn • e y = 0 


(A. 7) 
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Since t ^ 0 , from equation (A. 5), one concludes that 


Theorem 2 


Proof: 


= at 


(A.8) 


where a is a scalar function. Equations (A. 6) to (A— 8 ) imply that 
t • va + (v • t)a = 0 universally (A. 9) 


and 


o=0 at the inlet surface. 


(A. 10) 


From equations (A. 9) and (A. 10) and the assumption t ^ 0, we con- 
clude that a = 0 universally (ref. 17). As a result, 6a = at = 0 
universally. Thus o' = a" Q.E.D. 

: Assuming t ^ 0 and either vh ^ 0 or vS £ 0 throughout the 
region of interest, then for any solution n of equattions (A.l) 
and (A. 2), there exist infinitely many pairs of functions x and 
u (referred to as the Clebsch potentials associated with a) such 
that 

fl = vh Q x vt + vS x vp (A. 11) 


t • vx = 1 (A. 12) 

and 

t • v u = -T (A. 13) 


Conversely, tor any x and y which satisfy equations (A. 12) and 
(A. 13), equation (A. 11) defines a solution a of equations (A.l) 
and (A. 2). 

With the help of equations (A. 3), (A. 4), (A. 12), and (A. 13), the 
second part of this theorem can be proved by directly substituting 
equation (A. 11) into equations (A.l) and (A. 2). To prove the first 
part, first we assume that vh ^ 0. For any solution n of equa- 
tions (A.l) and (A. 2), one introduces a vector field o' such that 

a ' = a - vS x vp (A. 14) 

Here y is any function which satisfies equation (A. 13). Equa- 
tions (A.l), (A. 2), (A. 4), (A. 13), and (A. 14) imply that 
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(A. 15) 


^ x n' = vh Q 


and 


7 • ft' = 0 


(A. 16) 


As a result of equation (A. 15), one has 

n' • vh„ = 0 
o 


(A. 17) 


Equation (A. 17) states that h Q is conserved along any vortex line 

of o'. Using equations (A— 16 ) and (A. 17), it is shown in refer- 
ence 18 that there is a scalar function t such that 

o' = vh Q x 7 t (A. 18) 


provided that vh Q 4 0. Substituting equatiaons (A. 18) into (A. 15) 

and using equation (A. 3), one obtains equation (A. 12) as the con- 
dition t must satisfy. Equation (A. 11) follows directly from 
equations (A. 14) and (A. 18). Since there are infinitely many 
choices for y which satisfies equation (A. 13), it is easy to see 
that there are infinitely many pairs of Clebsch potentials asso- 
ciated with an n which satisfies equations (A.l) and (A. 2). 

If vh Q = 0, but vS 4 0, a similar proof can be constructed 

to show that the first part of Theorem 2 is still valid Q.E.D. 

The above theorem indicates that there^are infinitely many pairs of 
Clebsch potential associated with an n which satisfies equations 
(A.l) and (A. 2). The relationship among these Clebsch potentials 
is established by Theorem 3. 

Theorem 3: If t and y are a pair of Clebsch potentials associated with an 

n which satisfies equations (A.l) and (A. 2), then in a region in 
which vh Q x vS 4 0, t' and y' form another pair of Clebesn 

potentials associated with the same n if and only if there is a 
function q>(h Q ,S) such that 


3 CD ( h . S ) 



and 
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(A. 20) 


Proof : 


3<p(h ,S) 
-“ + — ? 5 — 


By directly substaituting equations (A. 19) and (A. 20) into equa- 
tions (A. 11) to (A. 13) and using equations (A. 3) and (A. 4), it can 
be shown that (t',u') and (t,u) are indeed associated with the same 
n if they are related by equations (A. 19) and (A. 20). To show the 
"only if" part of this theorem, first we define 

6t — t 1 - t and 6y = y‘ - y (A. 21) 


Since, by assumption, both (t',u') and (x,y) are associated with 
the same n, equations (A. 11) to (A. 13) and (A. 21) imply that 

7 h Q x 7 6 x + vS X 76 y = 0 (A. 22) 


K • 76 t = 0 (A. 23 ) 

and 

1 / • 76 y = 0 (A. 24 ) 

Using equation (A. 22), one has 

7 S • 7 h Q x 7 6 t = 0 (A. 25) 


and 


7 h Q • 7 S x 76 y = 0 


(A. 26) 


Equations (A. 25) and (A. 26) used in conjunction with the assumption 
that 7h Q x 7S ^ 0 imply that (ref. 19) 

6x = fi( h 0 >S) and 6y = f 2 (h Q ,S) (A. 27) 


Here fj and f^ are two arbitrary functions of h Q and S. 

Substituting equations |A.27)^into equation (A. 22) and again in- 
voking the assumption 7h Q x vS £ 0 yields 


3f^ 3f 2 

3 S = TFT 


(A. 28) 
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As a result of equation (A. 28), there exists a function <p(h ,S) 
tnat 0 

3q>(h ,S) . . 

= f 1 (* 1 0 »S) = — (A. 29) 


a<p(h ,S) 

5m = f 2 (h o ,S) = £ (A. 30) 

6t and 5y as given by equations (A. 29) and (A. 30) automatically 
satisfy equations (A. 23) and (A. 24) if equations (A. 3) and (A. 4) 
are taken into account. Equations (A. 19) and (A. 20) follow 
directjy from equations (A. 21), (A. 29), and (A. 30) Q.E.D. Since 
vh Q x vS = 0 for both homentropic and homenenergetic flows. 

Theorem 3 cannot be applied to these flows. Therefore the follow- 
ing theorems are given: 

Theorem 4: In a region in which vS = 0 and vh Q ^ 0, the function m intro- 

duced in Theorem 2 need not be specified (see eq; (A. 11)). If t 
is a Clebsch potential associated with an n which satisfies equa- 
tions (A.l) and (A.2),^then t* is another Clebsch potential asso- 
ciated with the same n if and only if there is a function y(h ) 
such that " 

x 1 = t + y( h Q ) (A. 31) 


Theorem 5: In a region in which vh Q = 0 and vS = 0, the function T intro- 
duced in Theorem 2 need not be specified.. If m is a Clebsch 
potential associated with an n which satisfies equations (A.l) 
and (A. 2), then m 1 is another Clebsch potential associated with 
the same ,n if and only if there is a function e(S) such that 

m 1 = m + 0 (S) (A. 32) 

» » » 4 

Theorems 4 and 5 can be proved using the same techniques employed 
in the Proof of Theorem 3. The .details will not be presented. 
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APPENDIX B 


Inlet Flow Conditions, Outflow Conditions, and Boundary Conditions 

In this appendix, we restrict the analysis to flows which enter the 
curved channel descrioed in Section III with parallel streamlines. Such flows 
have inlet conditions described by the following equations: 

V 1 = r(x 2 ,x 3 ), V 2 = V 3 = 0 (B.l) 


P = 


(B.2) 


and 


T = q(x 2 ,x 3 ) (B.3) 

2 3 

where P is a constant and r,q are arbitrary functions of x and x . 

(As noted in Section III, far upstream and downstream of the channel bend, the 
values of the metric tensor and thus the forms of fluid dynamic equations are 
identical in computational and physical space. As a result, equations (B.l) 
to (B.3) are valid in both spaces.) Using equations (B.l) to (B.3), it can be 
shown that, at the inlet 


T o ■ « + k 


(B.4) 


P o= P c( 1 + 2T 


Y 

Y-l 


(B.5) 


With the help of equations (II. 1), (II. 3), and (II. 4), equations (B.l), (B.4), 
and (8.5) can be used to establish the inlet conditions for the substitute 
t 1 ows . 

Since the streamlines of the incoming flow are parallel, the inlet stream- 
wise vorticity a 1 must vanish. From equation (II. 1), the inlet streamwise 

i 1 

vorticity a of the corresponding substitute flow must also vanish. Accord- 
ing to equation (1.19), a choice for the inlet profile of t' and y' which 
1 1 

yields a = 0 at the inlet is: 


t' = y' = 0 (B.6) 

With the inlet profiles of t' and y' established, the solution of their 
respective transport equation is uniquely defined. 

The construction of the appropriate inlet boundary conditions (B.C.S.) tor 

q>^, and is based on the assumption that the inlet flow condi- 
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tions do not vary during the course of iteration. As a result, the iterative 
increments f^ n ^, fj n ^, and ?j n ^ must satisfy the conditions: (n = 1 , 2 ,...) 

fi(n) = f(n) = j{n) _ Q (x 1 = a ) (B.7) 

A choice of the inlet B.C.S. for <|i^, and which is consistent 

with both the above condition and the gauge condition (11.32) is given in the 
in the first column of Table III. 

At the solid walls of the channel, we require the velocity component nor- 
mal to the walls to vanish. This requirement may be expressed in the tensor 
form 


V i n 1 = 0 


(B. 8 ) 


where and n 1 are, respectively, the covariant velocity vector and the 

contravariant normal vector at any point on the wall. For the numerical exam- 

2 2 

pie treated in Section III, the solid walls are located at x = c, x = d, 

3 3 

.x = 0, and x = e (See fig. 4). From equations (B. 8 ), (II. 10), and the fact 

that g. . = 0 if i 4 j, the boundary conditions at these walls become 

* J 

V 2 = F 2 = 0 at x 2 = c or x 2 = d (B.9) 

V 3 = F 3 = 0 at x 3 = 0 or x 3 = 0 (B.10) 


If we initialize the iteration procedure with a flow field which is consistent 
with the above conditions, the iterative increments f 1 ^, f^, and f^ 

(i = 2, 3) must satisfy the following conditions: (n = 1, 2, ...) 


f 2 (") „ f (n) . T (n) . 0 

f 3 ( n > = f . f ■ 0 


at 

at 




(B.ll) 

(B.12) 


A choice of the B.C.S. for 9 ^, \|»^, and at the solid walls which is 

consistent with both the above conditions and the gauge condition (11.32) is 
given in the second and third columns of Table III. 

At the channel exit, a fixed set of outflow conditions are directly im- 
posed on the substitute flow in computational space. The initial outflow con- 
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ditions are identical to the initial inlet flow conditions. The exit B 

for q/ n \ 4 >!j n \ an d are given in the fourth column of Table III 

B.C.S. insure the gauge condition equation (11.32) is satisfied, and in 
limit of incompressible flow yield 


for any n. 


C.S. 

These 

the 

(B.13) 
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APPENDIX C 


Inner Loop Stability Analysis 

In this appendix, the stability of the inner loop iteration is analyzed 
assuming the density p is held fixed. This assumption would be compatible 
with a solver in which the density is updated only in the outer loop. 

As in Section III, we shall use tensor notations in this study to econ- 
omize the mathematical expressions even though they may not be tensor expres- 
sions. We assume that v| n ^ and F 1 ^ are updated by equations (11.20) and 
(11.21) if n = 1, 3, 5, ... and updated by equations (11.22) and (11.23) if 
n = 2, 4, 6, .... Furthermore, we assume that F^ n ^ and are solved 

with the help of equations (11.27) and (11.33). (V^ could also be solved 

using equation (11.41) if n = 3, 5, 7, .... This alternative approach, how- 
ever, does not alter the resulting convergence criterion equation ( C . 11 ) ) . The 
above assumptions coupled with equations (11.12), (11.17) to (11.19), (11.28), 
and (11.34) imply that, for n = 2, 4, 6, ... 


and 


3\j; 


(n+1) 


3X° 3X 


J 


= e 


i jk 




3 y 


(n) 


ax J \ p Vg 3 x A 


ax 1 ax 1 



Vg g ij e jla 



(C.l) 


(C.2) 


To study the stability of the inner loop iterations, we shall assume that 


(n) ^ _(”) 


qr and ^ 


(n+1) ,(«■>) 


* 


as n -► +<= 


( C. 3) 


Furthermore, it is assumed that (n = 2, 4, 6, ...) 


and 





- Re(^ n )e IP J xJ ) 

Re(i' n+ 1 ) e IPjXJ ) 


(C.4) 


(C.5) 
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where P is the wave propagation vector, and ) the perturbation 

J *™’l 

amplitudes. Re = the real part of the quantity in the parenthesis, and 
I = 

This assumption is realistic only if the magnitude of is large (i.e., 
we deal with only high wave number disturbance) compared with the characteris- 
tic wave numbers of the spatial variations of p, g, g i ., and g 1J . As a re- 
suit, p, g d , and g. . are considered as constants in equations (C.l) and 

' J 

(C.2). A direct substitution of equations (C.3) to (C.5) into equations (C.l) 
and (C.2) yields 


and 


ii 


(n+l) 


i jk 


V^9 


g. P .p m ^ n ^ 
a k«, j jj- 


£ (" +2 > . vg 


(C.6) 


(C.7) 


Here P. are the directional cosines of the vector P., i.e., 

i i 


p def. _i 


fTi 


and thus 


■ ( ? i ) 2 * ( p 2 ) 2 + ( ? 3 ) 2 - 1 


(C.8) 


(C.9) 


Note that equations (C.5) and (C.6) are consistent with the gauge condition 
equation (11.32). If one substitutes equation (C.6) into equation (C.7), the 
final result is 


(n + 2) 

3L 


1 



(C.10) 


Since n = 2, 4, 6, ..., equation (C 

inner loop iterations requires that 
equivalently 


,10) implies that the stability of the 
X - < 1 ° r 
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2 > (g’Vj) (9JA)>° 


(C.ll) 


for any P which satisfies equation (C. 8 ). 

To solve equation (C.ll), we note that the matrices formed by giJ and 
9 ^ are positive definite hermitian matrices and are the inverse of each 
other. Assuming that x i, x?, and *3 Ul2 x 2 1 x 3 > are the eigenvalues 
of g^, it can be shown that (ref. 20 ): 



(C.12) 


/ (4XjX 3 ) and 1 are, respectively, the 
lower bound of ^g iJ P.,Pj^ (pki/k^s,) ^ the 
doman of includes all P which satisfy equation (C.9). Combining equa- 
tion (C.ll) with equation (C.12), one obtains 


O 

It should be noted that ^ + x^) 
least upper bound and the greatest 


6 



(C.13) 


For the current numerical study, g^ is given by equations (III. 5) and 
(III. 6 ). As a result, equation (C.13) reduces to 


2 

1_ > COS ( irX ) 1_ 

yn - ’ COSh(irX^) - v 2 


(C.14) 


The validity of the above stability condition has been verified by numerical 
experiments . 
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TABLE 1. - FLOW PROBLEMS AND THEIR DEFINING PARAMETERS 


Problem 


Geometric 

parameters 


Flow parameters 



a 

* 

c 

d 

6 

Y 

ma 

n 

D 


i 


NX1S 

No. l a 

-2 

a 

2 

a 

a 

45 

0. 

55 

0 

1 

i 

4 

0.1 

0.01 

0 

0 

72 

6 

6 


No. 2 a 













0.1 

0.02 



72 

6 

6 


1 

No. 3 b 













0.4733 

0.2 



144 

12 

12 


1 

No. 4 b 












, 

0.4733 

0.2 



108 

9 

9 

1 

1 

No. 5 b 













0.4733 

0.2 



72 

6 

6 

I 

I 

No. 6 C 













0.4733 

0.0 

0 

3 

144 

12 

12 

1 

1 


a Homentropic and linear. 

^Homentropic and nonlinear. 
c Inhomentropic, inhomeneryetic, and nonlinear. 


TABLE 2. - NUMERICAL COMPARISONS WITH EQUATION (111.16) 



Station 

A 

8 

C 

D 

E 

e (radian) 

0.6750 

1.5708 

2.4667 

2.8965 

3.0552 

Problem No. 1 

(n 0 = 0.1) 

5 

0.013432 

0.031329 

0.049137 

0.057498 

0.060317 

C/(2n 0 9) 

0.99498 

0.99722 

0.99604 

0.99254 

0.98712 

Problem No. 2 

(n 0 = 0.2) 

£ 

0.026856 

0.062547 

0.097549 

0.11294 

0.11708 

t/(2n 0 e) 

0.99470 

0.99546 

0.98870 

0.97475 

0.95806 





















































































If con verged, leave the outer loop 



Figure 1. - The flow chart for munk-prim flow solver. 























|b) Station C. 
Figure 6. - Concluded. 




(a) Station B. 

Fiqu re 7. - Secondary velocity field for problem (3). 



95) Station C. 


Figure 7. - Continued. 
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(c) Station E. 
Figure 7. - Concluded. 



(a) Station B. 

Figure 8. - Total pressure contours for problem (3). 




(b) Station C. 
Figure 8. -Continued. 



(c) Station E. 
Figure 8. - Concluded. 













Problem (3) 
Problem (4) 
Problem (5) 






to) Station C. 
Figure 10 - Continued. 



(c) Station E. 
Figure 10 - Concluded. 
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